%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                               Figures                                   %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% --------------------------- Description --------------------------------- 

    % This script plots the charts that are displayed in the main text of
    % the article.

%--------------------------------------------------------------------------
%                        Transitional Dynamics                            %
%--------------------------------------------------------------------------

    % Set the number of years before the shock
    pre_y = 30;
    add_pts = pre_y/p.yearlength - 1;

    % Population growth rate path
    eta_vector = InfoMatrix(1,:);
    eta_plot = [eta_vector(1)*ones(1,add_pts), eta_vector];

    % Population to varieties ratio path
    LFN_t = InfoMatrix(2,:);
    LFN_plot = [LFN_t(1)*ones(1,add_pts), LFN_t];

    % Variety growth rate path
    gN_vector = InfoMatrix(3,:);
    gN_plot = [gN_vector(1)*ones(1,add_pts), gN_vector];

    % Production share path
    prod_share_vector = InfoMatrix(4,:);
    prod_share_plot = [prod_share_vector(1)*ones(1,add_pts), prod_share_vector];

    % Misallocation Term path
    M_vector = InfoMatrix(5,:);
    M_plot = [M_vector(1)*ones(1,add_pts), M_vector];

    % Computes the growth rates of the Misallocation term
    gMis = [0, log(M_vector(2:end)./M_vector(1:end-1))];
    gMis = movmean(gMis, [(1/p.yearlength-1) 0]) .* min(linspace(1, TimePeriods, TimePeriods), 1/p.yearlength);
    gMis_plot = [gMis(1)*ones(1,add_pts), gMis];

    % Labor Wedge Term path
    Lam_vector = InfoMatrix(6,:);
    Lam_plot = [Lam_vector(1)*ones(1,add_pts), Lam_vector];    
   
    % Compute the productivity growth rate along the transition
    gQ_plot = -(p.qbar - 1)*(gN_plot + p.d0)/(1-p.alpha);
    g_plot = -gQ_plot/(p.s-1) + p.In + gN_plot /(p.s-1) + gMis_plot;

    % Firm moments equilibrium paths
    Size_plot = [AvFirmSize_t(1)*ones(1,add_pts), AvFirmSize_t];
    Entry_plot = [EntryRate_t(1)*ones(1,add_pts), EntryRate_t];
    Exit_plot = [ExitRate_t(1)*ones(1,add_pts), ExitRate_t];

    % U function equilibrium path
    ulam_plot = u_func(1,:);
    ulam_plot = [ulam_plot(1)*ones(1,add_pts), ulam_plot];
    uces_plot = u_func(end,:);
    uces_plot = [uces_plot(1)*ones(1,add_pts), uces_plot];    

    % Set the time grid to the year standard
    t_vector = p.TimeGrid*(p.Years/p.TimeGrid(end)) + 1980;
    t_plot = [t_vector(2:pre_y/p.yearlength) - pre_y, t_vector];

    % ----------------- Fixed Point check ---------------------------------

    % Production share implied by market clearing condition
    prod_share_MC = 1 - ((gN_vector + p.d0)/(1-p.alpha) - (p.zeta-1)/p.zeta*p.xstar)./(p.phie.*LFN_t);
    prod_share_MC = [prod_share_MC(1)*ones(1,add_pts), prod_share_MC];

    % Derivative of u function wrt time
    ulambda = u_func(1,:);
    ulambda_d = [0, diff(u_func(1,:))*p.yearlength];
    uCES = u_func(end,:);
    uCES_d = [0, diff(u_func(1,:))*p.yearlength];

    % Production share growth rate
    gShare = [0, log(prod_share_vector(3:end)./prod_share_vector(2:end-1)), 0];
    gShare = movmean(gShare, [0 (1/p.yearlength-1)]) .* fliplr(min(linspace(1, TimePeriods, TimePeriods), 1/p.yearlength));
    gShare(1) = 0;
    
    % Recovers the growth of Lambda from InfoMatrix
    gLam = InfoMatrix(7,:);

    % Compute the term associated with u(lambda) and u(mu)
    lambda_term = p.alpha*p.l^(p.s-1)*(ulambda.*(p.rho-gLam-eta_vector+(gN_vector+p.d0)/(1-p.alpha))*p.yearlength - ulambda_d).*p.phie.*LFN_t;
    CES_term = (1-p.alpha)*p.wbar*(uCES.*(p.rho-gLam-eta_vector+(gN_vector+p.d0)/(1-p.alpha))*p.yearlength - uCES_d).*p.phie.*LFN_t;

    % Solve for the resulting production share in period t
    prod_share_free = (p.rho + gShare - gLam + (p.alpha*gN_vector+p.d0)/(1-p.alpha)  - (p.zeta-1)/p.zeta*p.xstar)./(lambda_term+CES_term);
    prod_share_free = [prod_share_free(1)*ones(1,add_pts), prod_share_free];

%--------------------------------------------------------------------------
%                           Control Charts                                %
%--------------------------------------------------------------------------

    % Production share path
    figure('Name','Transitional Dynamics - Control Figures');
    subplot(2,3,1)
    plot(t_plot, prod_share_plot*100, 'DisplayName','Model', 'LineWidth', 2)
    hold on
    yline(prod_share_counter*100,'--','DisplayName','BGP', 'Color', 'k');
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Production Share (%)", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Production Share", "FontSize", 14)
    legend('Location', 'southeast')

    subplot(2,3,2)
    plot(t_plot, prod_share_MC*100, 'linewidth', 1,'DisplayName','Mkt. Clearing');
    hold on
    plot(t_plot, prod_share_free*100, 'linewidth', 2, 'LineStyle', '--', 'DisplayName','Free Entry');
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Production Share (%)", "Rotation", 90, 'FontWeight', 'bold')
    title("Fixed Pt. Check", "FontSize", 14)
    legend('Location', 'southeast')

    subplot(2,3,3)
    plot(t_plot, ulam_plot*p.yearlength, 'DisplayName','Model', 'LineWidth', 2)
    hold on
    yline(k_func_lam_counter/(M_counter^(p.s-1)*Lam_counter^p.s),'--','DisplayName','BGP', 'Color', 'k');
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("b_t(\lambda)", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("b_t(\lambda)", "FontSize", 14)
    legend('Location', 'southeast')

    subplot(2,3,4)
    plot(t_plot, uces_plot*p.yearlength, 'DisplayName','Model', 'LineWidth', 2)
    hold on
    yline(k_func_mu_counter/(M_counter^(p.s-1)*Lam_counter^p.s),'--','DisplayName','BGP', 'Color', 'k');
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("b_t(\mu)", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("b_t(\mu)", "FontSize", 14)
    legend('Location', 'southeast')

    subplot(2,3,5)
    plot(t_plot, M_plot, 'DisplayName','Model', 'LineWidth', 2)
    hold on
    yline(M_counter,'--','DisplayName','BGP', 'Color', 'k');
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Misallocation Term", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Misallocation Term", "FontSize", 14)
    legend('Location', 'northeast')

    subplot(2,3,6)
    plot(t_plot, Lam_plot, 'DisplayName','Model', 'LineWidth', 2)
    hold on
    yline(Lam_counter,'--','DisplayName','BGP', 'Color', 'k');
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Labor Wedge Term", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Labor Wedge", "FontSize", 14)
    legend('Location', 'northeast')

%--------------------------------------------------------------------------
%                               Figure 12                                  %
%--------------------------------------------------------------------------

    % Computes the productivity term associated to quality growth
    gQ_term = -gQ_plot/(p.s-1) + p.In;

    % Computes the productivity term associated to variety growth rate
    gN_term = gN_plot /(p.s-1);

    % Productivity growth rate breakdown
    figure('Name','Productivity Growth Rate Breakdown');
    plot(t_plot, gQ_term*100, 'DisplayName', 'Efficiency Growth Term', 'LineWidth', 2)
    hold on
    plot(t_plot, gN_term*100, 'DisplayName', 'Variety Growth Term', 'LineWidth', 2)
    hold on
    plot(t_plot, g_plot*100, 'DisplayName', 'Productivity Growth', 'LineWidth', 2)
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Growth rates (%)", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Productivity Growth Breakdown", "FontSize", 14)
    legend('Location', 'northeast')

%--------------------------------------------------------------------------
%                               Figure 11                                  %
%--------------------------------------------------------------------------

    % Productivity growth rate
    figure('Name','Transitional Dynamics');
    subplot(1,2,1)
    plot(t_plot, g_plot*100, 'DisplayName','Model', 'LineWidth', 2)
    hold on
    yline(g_counter*100,'--','DisplayName','BGP', 'Color', 'k');
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Productivity growth rate (%)", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Productivity Growth", "FontSize", 14)
    legend('Location', 'northeast')

    % Variety Intensity ratio
    subplot(1,2,2)
    plot(t_plot, 1./LFN_plot, 'DisplayName','Model', 'LineWidth', 2)
    hold on
    yline(1/LFN_counter,'--','DisplayName','BGP', 'Color', 'k');
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Variety intensity", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Variety Intensity", "FontSize", 14)
    legend('Location', 'southeast')

%--------------------------------------------------------------------------
%                              Figure 9/10                                 %
%--------------------------------------------------------------------------

    plot10y = [Firms_10y_t(1)*ones(1,add_pts), Firms_10y_t];
    Prods2_plot = [Prods2(1)*ones(1,add_pts), Prods2];
    Prods5_plot = [Prods5(1)*ones(1,add_pts), Prods5];
    avgmu_plot = [AvgMup(1)*ones(1,add_pts), AvgMup];

    % Share of firms above 10 years
    figure('Name','Transitional Dynamics');
    subplot(2,2,1)
    plot(t_plot, plot10y*100, 'DisplayName','Model', 'LineWidth', 2);
    hold on
    yline(Firms_10y_counter*100, '--')
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Share of Firms above 10 years (%)", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Shifting Age Distribution", "FontSize", 14)
   
    % Number of products per firm
    subplot(2,2,2)
    yyaxis left
    plot(t_plot, Prods2_plot*100, 'DisplayName','2+ Products', 'LineWidth', 2)
    hold on
    ylabel("Share of Firms (%)", "Rotation", 90, 'FontWeight', 'bold')
    yyaxis right
    plot(t_plot, Prods5_plot*100, 'DisplayName','5+ Products', 'LineWidth', 2)
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Share of Firms (%)", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Rising Product Concentration", "FontSize", 14)
    legend('Location', 'southeast')

    % Unconditional expected markup
    subplot(2,2,3)
    plot(t_plot, avgmu_plot, 'DisplayName','2% Population growth', 'LineWidth', 2)
    hold on
    yline(CostMarkup_counter, '--')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Average Markup", "Rotation", 90, 'FontWeight', 'bold')
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Average Markup", "FontSize", 14)

    % Expected markup by firm age - two BGPs
    subplot(2,2,4)
    plot(avec, CostMarkup_A, 'DisplayName','2% Population growth', 'LineWidth', 2)
    hold on
    plot(avec, CostMarkup_A_counter, 'DisplayName','0.24% Population growth', 'LineWidth', 2);
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Average Markup", "Rotation", 90, 'FontWeight', 'bold')
    ylim([1 1.4])
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Markups Across Firms", "FontSize", 14)
    legend('Location', 'southeast')

%--------------------------------------------------------------------------
%                               Figure 8                                  %
%--------------------------------------------------------------------------

    % Population growth rate and variety growth rate
    figure('Name','Transitional Dynamics');
    subplot(2,2,1)
    plot(t_plot, gN_plot*100, 'DisplayName','Variety', 'LineWidth', 2)
    hold on
    plot(t_plot, eta_plot*100, 'DisplayName','Population', 'LineWidth', 1.5)
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Growth rates (%)", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Population and Variety Growth", "FontSize", 14)
    legend('Location', 'northeast')

    % Entry rate 
    subplot(2,2,2)
    plot(t_plot, Entry_plot*100, 'DisplayName','Model', 'LineWidth', 2)
    hold on
    yline(EntryRate_counter*100,'--','DisplayName','BGP', 'Color', 'k');
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Entry rate (%)", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Entry Rate", "FontSize", 14)
    legend('Location', 'northeast')

    % Exit rate 
    subplot(2,2,3)
    plot(t_plot, Exit_plot*100, 'DisplayName','Model', 'LineWidth', 2)
    hold on
    yline(ExitRate_counter*100,'--','DisplayName','BGP', 'Color', 'k');
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Exit rate (%)", "Rotation", 90, 'FontWeight', 'bold')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Exit Rate", "FontSize", 14)
    legend('Location', 'northeast')

    % Average Firm Size 
    subplot(2,2,4)
    plot(t_plot, Size_plot, 'DisplayName','Model', 'LineWidth', 2)
    hold on
    yline(AvFirmSize_counter,'--','DisplayName','BGP', 'Color', 'k');
    hold on
    xline(1980, '--', 'HandleVisibility', 'off', 'Color', 'r')
    xlabel("Year", "Rotation", 0, 'FontWeight', 'bold')
    ylabel({'{\bf Average Firm Size}'; '{Employees}'}, 'FontWeight', 'normal')
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    title("Average Firm Size", "FontSize", 14)
    legend('Location', 'southeast')

 %--------------------------------------------------------------------------
%                           Figures 6 and 7                               %
%--------------------------------------------------------------------------

    % Determines an age cut for the plots
    c = 240;

    % Reestablish the age vector used for these plots
    avec_plot = linspace(0, p.maxTime, p.maxTime/p.periodlength);

    % This chart represents Sales Growth by Age
    figure('Name','Calibration');
    subplot(2,2,1)
    plot(avec_plot(1:c), Esales(1:c), 'linewidth', 2, 'DisplayName','Calibrated Model');
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    xlabel("Firm Age", "Rotation", 0, 'FontWeight', 'bold')
    ylabel({'{\bf Log Average Sales}'; '{Difference from age 0}'}, 'FontWeight', 'normal')
    title("Sales Growth by Age", "FontSize", 14)
    legend('Location', 'southeast')

    % This chart represents Markup Growth by Age
    markuplot = Emarkup - Emarkup(1);
    subplot(2,2,2)
    plot(avec_plot(1:c), markuplot(1:c)*100, 'linewidth', 2, 'DisplayName','Calibrated Model');
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    xlabel("Firm Age", "Rotation", 0, 'FontWeight', 'bold')
    ylabel({'{\bf Markup}'; '{Pct. Points Above Age 0}'}, 'FontWeight', 'normal')
    title("Markup Growth by Age", "FontSize", 14)
    legend('Location', 'southeast')


    % This chart plots the employment distribution in discrete bins
    subplot(2,2,3)
    X = {'a) 1 to 4','b) 5 to 9','c) 10 to 19','d) 20 to 49','e) 50 to 99','f) 100 to 249','g) 250 to 499','h) 500 to 999','i) 1,000 to 2,499','j) 2,500 to 4,999','k) 5,000 to 9,999','l) 10,000+'};
    bar((categorical(X)), [empshare']*100);
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    ylabel("Employment Share (%)", "Rotation", 90, 'FontWeight', 'bold')
    title("Employment Distribution", "FontSize", 14)
    legend({'Calibrated Model'}, 'Location','northwest')

    % This chart plots the exit rate by age
    subplot(2,2,4)
    plot(avec_plot(1:c), Exit_age(1:c)*100, 'linewidth', 2, 'DisplayName','Calibrated Model');
    ax = gca;
    ax.XGrid = 'off';
    ax.YGrid = 'on';
    xlabel("Firm Age", "Rotation", 0, 'FontWeight', 'bold')
    ylabel("Lifecycle Exit Rate (%)", "Rotation", 90, 'FontWeight', 'bold')
    title("Exit Rates by Age", "FontSize", 14)
    legend

    % Clear the variables used to construct the charts
    clear X ax Emarkup_orig Esales_orig empshare_orig Exit_age_orig paramgues_orig markuplot markuplot_orig
    %clear add_pts avec_plot dvdt Entry_plot eta_plot etaseq Exit_plot g_plot gN_plot gQ_plot gShare i LFN_plot plot10y prod_share_plot
    %clear Prods2_plot Prods5_plot Size_plot t_plot ufun_d ulam_plot u_term prod_share_MC prod_share_free pre_y c

